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Abstract. 

The diffusivity of tagged particles is demonstrated to be heterogeneous on time 
scales comparable to or less than the structural relaxation time in a highly supercooled 
liquid via 3D molecular dynamics simulation. The particle motions in the relatively 
active regions dominantly contribute to the mean square displacement, giving rise to 
a diffusion constant systematically larger than the Stokes-Einstein value. The van 
Hove self-correlation function G s (r, t) is shown to have a large r tail which can be 
scaled in terms of r/t 1 / 2 for t < 3t q , where r a = the stress relaxation time. Its 
presence indicates heterogeneous diffusion in the active regions. However, the diffusion 
process eventually becomes homogeneous on time scales longer than the life time of 
the heterogeneity structure (~ 3r a ). 



INTRODUCTION 

Molecular dynamics (MD) simulations can be powerful tools to gain insights into 
relevant physical processes in highly supercooled liquids. In particular, a number 
of recent MD simulations have detected dynamic heterogeneities in supercooled 
model binary mixtures [1-6]. That is, rearrangements of particle configurations 
in glassy states are cooperative, involving many molecules, owing to configuration 
restrictions. Recently, we succeeded in quantitatively characterizing the dynamic 
heterogeneities in two (2D) and three dimensional (3D) model fluids via MD simu- 
lations. We examined bond breakage processes among adjacent particle pairs and 
found that the broken bonds in an appropriate time interval (~ the stress relax- 
ation time or the structural a relaxation time r a ) are very analogous to the critical 
fluctuations in Ising spin systems with their structure factor being excellently fit- 
ted to the Ornstein-Zernike form [3,4]. The correlation length £ thus obtained is 
related to r a via the dynamic scaling law, r a ~ £ z , with z = 4 in 2D and z = 2 
in 3D. The heterogeneity structure in the bond breakage is essentially the same as 
that in jump motions of particles from cages or that in the local diffusivity, as will 
be discussed below. In this paper, we investigate heterogeneities of tagged particle 
motions in a 3D supercooled liquid [5]. 

In a wide range of liquid states, the Stokes-Einstein relation Drja/ksT = const. 
has been successfully applied between the translational diffusion constant D of a 



tagged particle and the viscosity r] even when the tagged particle diameter a is of 
the same order as that of solvent molecules. However, this relation is systematically 
violated in fragile supercooled liquids [4,5,7-11]. The diffusion process in super- 
cooled liquids is thus not well understood. In particular, Sillescu et al. observed 
the power law behavior D oc r\~ v with v = 0.75 at low temperatures [8]. Further- 
more, Ediger et al. found that smaller probe particles exhibit a more pronounced 
increase of Drj/T oc D/Dse with decreasing T [9], where D$e ~ ksT /2irria is the 
Stokes-Einstein diffusion constant. In such experiments the viscosity changes over 
12 decades with decreasing T, while the ratio D/D SE increases from order 1 up 
to order 10 2 ~ 10 3 . The same tendency has been detected by molecular dynam- 
ics simulations in a 3D binary mixture with N = 500 particles [10] and in a 2D 
binary mixture with N = 1024 [11]. In our recent 3D simulation with N = 10 4 
[4,5], i] and D both varied over 4 decades and the power law behavior D oc ^~ - 75 
has been observed. Many authors have attributed the origin of the breakdown to 
heterogeneous coexistence of relatively active and inactive regions, among which 
the local diffusion constant is expected to vary significantly [8,9,12-14]. The aim of 
this paper is to demonstrate via MD simulation that the diffusivity of the particles 
is indeed heterogeneous on time scales shorter than r a but becomes homogeneous 
on time scales much longer than r a . 

SIMULATION METHOD AND RESULTS 

Our 3D binary mixture is composed of two atomic species, 1 and 2, with N± = 
N 2 = 5000 particles with the system linear dimension L = V 1 ^ 3 being fixed at 
23. 1o\ [15]. They interact via the soft-core potentials v a b(r) = e(a a f ) /r) 12 with 
Cab = (°a + 0&)/2, where r is the distance between two particles and a,b G {1,2}. 
The interaction is truncated at r = 3<7i. The mass ratio is m 2 /m 1 = 2. The 
size ratio is o 2 ja\ = 1.2, which is known to prevent crystallization [1,15]. No 
tendency of phase separation is detected at least in our computation times. We fix 
the particle density at a very high value of (Ni + N 2 )/V = 0.8/a 3 , so the particle 
configurations are severely restricted or jammed. We will measure space and time 
in units of a\ and r = (mi of/ 'e) 1 / 2 . The temperature T will be measured in units 
of e/ks, and the viscosity 77 in units of er /af. The time step At = 0.005 is used. In 
our systems the structural relaxation time becomes very long at low temperatures. 
Therefore, very long annealing times (2.5 x 10 5 for T = 0.267) are chosen in our 
case. For T > 0.267, no appreciable aging (slow equilibration) effect is detected 
in various quantities such as the pressure or the density time correlation function, 
whereas at T = 0.234, a small aging effect remains in the density time correlation 
function. 

Let us first consider the incoherent density correlation function, F s (q,t) = 
(Ylfli exp [iq ■ Arj(t)])/Ni for the particle species 1, where Ar 7 (t) = rj(t) — rj(0) 
is the displacement vector of the j'-th particle. The a relaxation time r a is 
then defined by F s (q,r a ) = e _1 at q = 2n for various T. We also calcu- 
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FIGURE 1. Coherent and incoherent intermediate scattering functions for various temperatures 
with q = 2ir (a peak wavenumbcr in S n (q,0)). T = 0.234, 0.267, 0.306, 0.352, 0.473, and 0.772 
from right. 

late the coherent time correlation function, Sn(q, t) = {ni(q,t)ni(—q,0)), where 
ni(q, t) = J2f=i exp [iq ■ Vj(t)] is the Fourier component of the density fluctuations 
of the particle species 1. The decay profiles of Sn(q, t) at its first peak wave number 
q = Qm ~ 27T and F s (q, t) at q = 2tt nearly coincide in the whole time region studied 
(t < 2 x 10 5 ) within 5% as shown in Fig. 1. Hence Su(q m , T a )/Su(q m , 0) = e _1 
holds for any T in our simulation. Such agreement is not obtained for other wave 
numbers, however. These results are consistent with those for a Lennard- Jones bi- 
nary mixture [16]. Furthermore, some neutron-spin-echo experiments [17] showed 
that the decay time of Su(q m , t) is nearly equal to the stress relaxation time and as 
a result the viscosity r\ is of order r a . In agreement with this experimental result, 
we obtain a simple linear relation in our simulations [5], 

t« = (AJql^/T. (1) 

The coefficient A, q is close to 2-7T in our system. Here, we may define a (/-dependent 
relaxation time r q by F(q,T q ) = e -1 . Thus, particularly at the peak wave number 
q = q m , the effective diffusion constant defined by D q = l/q 2 r q is given by the 
Stokes-Einstein form even in highly supercooled liquids. However, notice that the 
usual diffusion constant is the long wavelength limit, D = \im q ^ D q . It is usually 
calculated from the mean square displacement, ((Ar(t)) 2 ) = (J2f=i(^ 3 (t)) 2 ) /N v 
The crossover of this quantity from the plateau behavior arising from motions in 
transient cages to the diffusion behavior 6Dt has been found to take place around 
t ~ 0.lT a [4]. In Fig.2, we plot Dr a versus r Q , which clearly indicates breakdown 
of the Stokes-Einstein relation in agreement with the experimental trend. 

To examine the diffusion process in more detail, we consider the van Hove self- 
correlation function, G s (r,t) = (Ef=i S{Arj(t) - r))/N v Then, 
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FIGURE 2. DT a versus r Q in a quiescent supercooled liquid. The solid line represents the 
Stokes-Einstein value DsE T a — (2tt)~ 2 arising from Eq.(l). 



r 

F s {q,t)= / 
jo 



J sm(gr) 2 . 



(2) 



is the 3D Fourier transformation of G s (r,t). At the peak wavenumber q = 2tt, 
the integrand in Eq.(2) vanishes at r = 1, and the integral in the region r < 1 is 
confirmed to dominantly determine the decay of F s (2ir,t). On the other hand, the 
mean square displacement 
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/ drAiir A G s {r,t) 
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(3) 



is determined by the particle motions out of the cages for t > r Q in glassy states. In 
Fig. 3, we display 47rr 4 G s (r, r a ) versus r, where r a = 3.2 and 2000 for T = 0.473 and 
0.267, respectively. These curves may be compared with the Gaussian (Brownian 
motion) result, (2/vr) 1 / 2 £- 3 r 4 exp(-r 2 /2£ 2 ), where U 2 = 6D SE T a = 3/2vr 2 is the 
Stokes-Einstein mean square displacement. Because the areas below the curves 
give QDr a , we recognize that the particle motions over large distances r > 1 are 
much enhanced at low T, leading to the violation of the Stokes-Einstein relation. 

We then visualize the heterogeneity of the diffusivity. To this end, we pick up 
mobile particles of the species 1 with |Arj(t))| > ^ c in a time interval [to,to + t] and 



number them as j 



N m . Here £ c is defined such that the sum of Arj(ty 



of the mobile particles is 66% of the total sum 
these particles are drawn as spheres with radii 



6DtNi for t > 0.1r Q ). In Fig.4, 



a s (t) = \Ar s (t)\/y/((Ar(t))*) 



(4) 



located at Rj(t) = h[fj(to) + r j(to + t)] i n time intervals [to, t + r a ] for T = 0.473 
(a) and 0.267 (b) and in [t ,t + 10r a ] for T = 0.267 (c). The the mobile 
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FIGURE 3. 47rr 4 G ;j (r, t) versus r at t = r a . The solid line is for T = 0.267 and the broken line 
is for T — 0.473. The dotted line represents the Brownian motion result. The peaks at r ~ 1.2 
and 2 in the solid line arise from hopping processes in our system at T = 0.267. Note that the 
areas below the curves give 6Dr a . 

particle number N m is 1595 in (a), 725 in (b), and 1316 in (c), respectively. 
Here the Gaussian results is N m = 1800. The ratio of the second moments 
c 2 = Ejh a i(*)7E^i%(*) 2 is held fixed at 0.66, while the ratio of the fourth 
moments c 4 = J2f=i a j(t) 4 / Hf=i CLj(t) 4 turns out to be close to 1 as c 4 = 0.89 
in (a), 0.92 in (b), and 0.90 in (c). The mobile particles are homogeneously 
distributed for T = 0.473 at r a , whereas for T = 0.267, the heterogeneity is 
significant at r a , but is much decreased at 10r Q . In fact, the variance defined 

by V = A^ m E^i%(0 4 /(E^iaj(0 2 ) 2 - 1 is °- 27 in ( a ), °- 41 in ( b ), and °- 32 
in (c). Note that the statistical average of V (taken over many initial times 
to) is related to the non-Gaussian parameter A 2 = 3(Ar(t) 4 )/5(Ar(t) 2 ) 2 — 1 = 
3iV 1 (Efii%(t) 4 )/5((Efiia j (t) 2 )) 2 -lby 

(V) ^ (5(c 4 )(iV m )/3c 2 iV 1 )(l + A 2 (t))-l, (5) 

where the deviations c 4 — (c 4 ) and N m / (N m ) — 1 are confirmed to be very small 
for large N\ and are thus neglected. We may also conclude that the significant 
rise of A%{t) in glassy states originates from the heterogeneity in accord with some 
experimental interpretations [18]. 

We next consider the Fourier component of the diffusivity density defined by 

Nm 

V q (t ,t) = J2 a j(t) 2 eM-i<l-RAt)], (6) 

which depends on the initial time to and the final time to + t. The correlation 
function Sv(q, t, r) = (Vq(t + r, t)P_q(t , t)) is then obtained after averaging over 





FIGURE 4. Mobile particles of the species 1 with a time interval t = T a at T = 0.473 
(a) and T = 0.267 (b) and with t = 10r Q at T = 0.267 (c). The radii of the spheres are 
\Arj(t)\/^ ((Ar(t)) 2 ) and the centers are at |[^-(to) + rj(t + t)]. The system linear dimension 
is L = 23.2. The darkness of the spheres represents the depth in the 3D space. 
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FIGURE 5. The correlation functions Sv(q, T a , 0) (a) and Sb{q, T a , 0) (b) with q = 2ir. 



many initial states. We plot Sx>(g, r a , 0) in Fig. 5 (a). The heterogeneity structure 
S&(9> T a , 0) of the bond breakage [3,4] with a time interval of r a is also plotted in 
Fig. 5 (b). It is confirmed that Sx>(g, r a ,0) tends to its long wavelength limit for 
q ^ where £ coincides with the correlation length obtained from Sf,(q, r a , 0). 
As the difference r of the initial times increases with fixed t = r a , Sx>{q,T a ,T) 
relaxes as exp[— {r/T h ) c } for q < where c ~ 0.5 at T = 0.267 and r h ~ 3r a 
is the life time of the heterogeneity structure. The two-time correlation function 
among the broken bond density [3,4] also relaxes with in the same manner. 

We naturally expect that the distribution of the particle displacement Arj (t) in 
the active regions should be characterized by the local diffusion constant D(x,t) 
dependent on the spatial position x = (x, y, z) and the time interval t. The van 
Hove correlation function G s (r,t) may then be expressed as the spatial average of 




10t q < t < IOOtq, (t = 10nT Q , n = 1,2, •••,10), where the heterogeneity effect is smoothed 
out. The dotted lines are the Gaussian form. 

a local function G s (x,r,t), which is given by [4irD(x, t)t}~ 3 ^ 2 exp[— r 2 /AD(x, t)t\. 
To check this conjecture, we plot the scaled function y 6Dt47rr 2 G s (r, t) versus r* = 
r/VODi in Fig. 6. The areas below the curves are fixed at 1. At relatively short 
times t ^ 3r a , the curves in the region r > 1 or r * > (6Dt)- 1/2 , which give dominant 
contributions to ((Ar(t)) 2 ), tend to a master curve quite different from the rapidly 
decaying Gaussian tail. Note that in each curve the position of the peak at larger 
r* corresponds to r = 1. This asymptotic law is consistent with the picture of the 
space-dependent diffusion constant in the active regions. It is also important that 
the heterogeneity structure remains unchanged in the time region t ^ ~ 3r Q . 
At longer times t > 10r a , the curves approach the Gaussian form as can be seen 
in the inset of Fig. 5. Of course, 47rr 2 G s (r,t) for r < 1 does not scale in the above 
manner, because it is the probability density of a tagged particle staying within a 
cage. This short-range behavior determines the decay of F s (2ir,t) as noted below 
Eq.(2). 

CONCLUDING REMARKS 

In our previous studies [3,4], we performed extensive MD simulations and identi- 
fied weakly bonded or relatively active regions from breakage of appropriately defined 
bonds. We also found that the spatial distributions of such regions resemble the 
critical fluctuations in Ising spin systems, so the correlation length £ can be deter- 
mined. It grows up to the system size as T is lowered, but no divergence seems to 
exist at nonzero temperatures. In the present work, we have demonstrated that the 
diffusivity in supercooled liquids is spatially heterogeneous on time scales shorter 



than 3r a , which leads to the breakdown of the Stokes-Einstein relation [5]. The 
heterogeneity detected is essentially the same as that of the bond breakage in our 
previous works [3,4]. We should then investigate how the heterogeneity arises and 
influences observable quantities in more realistic glass-forming fluids with complex 
structures. 
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